RNA_Sequencing

Loading the required libraries

library(tidyverse)      # for data manipulation
library(DESeq2)         # for Differential gene expression analysis
library(ggplot2)        # for plotting
library(plotly)         # for interactive plots
library(GenomicRanges)  # for creating and manipulating genomic ranges
library(mixOmics)       # for PCA

Importing data

assay<- read.csv("assay.csv", row.names = 1)  # Importing gene expression assay data

Warning in file(file, "rt") :
  cannot open file 'assay.csv': No such file or directory
Error in file(file, "rt") : cannot open the connection

Printing all data

Visualizing variability in library size

ggplotly(lsize)

Warning: The following aesthetics were dropped during statistical transformation: label.
ℹ This can happen when ggplot fails to infer the correct grouping structure in the data.
ℹ Did you forget to specify a `group` aesthetic or to convert a numerical variable into a factor?

Filtering out low count genes

Filtering low-count genes is a common preprocessing step in RNA-Seq data analysis, especially when using tools like DESeq2 for differential expression analysis. Here’s why it’s important:

Low-Count Genes Provide Little Statistical Power

Genes with low counts typically have high variability due to the limited sampling of sequencing reads. This variability makes it difficult to reliably detect significant differences in expression levels between conditions.

Reduction of Multiple Testing Burden

Differential expression analysis involves statistical tests for thousands of genes. Each test adds to the multiple-testing burden, requiring more stringent p-value adjustments.

Minimizing Noise

Low-count genes are often affected by: Sequencing noise: Random fluctuations during library preparation or sequencing. Mapping artifacts: Reads mapped spuriously to regions with minimal biological relevance.

We are usnig the following filters

rowSums(Exp.data > 0) >= 2: At least 2 samples with non-zero expression.

rowSums(Exp.data) >= 20: Total expression across samples is at least 20

paste0("Number of genes after filtering: ", nrow(assay.filtered))

[1] "Number of genes after filtering: 201"

Visualizing RLE

A Relative Log Expression (RLE) plot is a diagnostic tool commonly used in RNA-Seq data analysis to assess the normalization of a count matrix. It shows the log ratios of gene expression counts relative to their median across samples.

Creating DEseq2 object with filtered counts

In DESeq2, the primary data structure for storing RNA-Seq data and performing analyses is the DESeqDataSet (DDS) object. It is specifically designed to handle the raw counts and associated metadata for differential expression analysis.

What is a DESeqDataSet?

A DESeqDataSet is an object from the S4 class system in R that holds:

Count Data: The raw, unnormalized read counts for genes across samples.

Column Metadata: Information about the experimental design (e.g., condition, batch).

Row Metadata: Gene annotations or other feature-level metadata.

Design Formula: Specifies the experimental design to model differential expression (e.g., ~ condition or ~ condition + batch).

Note: DESeq2 take raw counts as input

DE.obj <- DESeqDataSetFromMatrix(countData = assay.filtered,
                                 colData   = colData,
                                 metadata  = rowRanges,
                                 design    = ~ Sample_Group)

Error in getListElement(x, i, ...) : 
  GRanges objects don't support [[, as.list(), lapply(), or unlist() at the moment

Estimating size factor and normalizing counts

DESeq2 performs normalization using ‘size factor’ to account for differences in sequencing depth and other systematic biases across samples in RNA-Seq datasets.

The size factor is a scaling factor used in RNA-Seq analysis to normalize read counts across samples, accounting for differences in sequencing depth or library size. This ensures that comparisons of gene expression levels are not biased by variations in the total number of reads per sample.

Visualizing data

Principal Component Analysis (PCA) is a dimensionality reduction technique widely used in data analysis, including RNA-Seq or other high-dimensional datasets. It transforms the data into a new coordinate system such that the greatest variance lies along the first principal component (PC1), the second greatest variance along the second component (PC2), and so on.

Performing Differential expression analysis

DESeq2 performs differential expression gene (DEG) analysis using a model-based approach designed specifically for count data from RNA-Seq experiments. It employs generalized linear models (GLMs) and statistical techniques to identify genes with significant differences in expression across conditions.

Normalization

Dispersion Estimation

Model fitting

Hypothesis Testing

Multiple Testing Correction

DE.obj$Sample_Group <- relevel(DE.obj$Sample_Group, ref = "HK")  # Setting reference group

DE.obj <- DESeq(DE.obj)                                          #Fitting the model

using pre-existing size factors
estimating dispersions
found already estimated dispersions, replacing these
gene-wise dispersion estimates
mean-dispersion relationship
final dispersion estimates
fitting model and testing

res.con <- results(DE.obj,                                       # Extracting result table
                   contrast=c("Sample_Group", "WT", "HK"),       # Setting Contrast 
                   pAdjustMethod = "BH",                         # Using 'BH' for p-value adjustment
                   alpha = 0.1)                                  # Adjusted p-value threshold
                   
res.con <- res.con[!is.na(res.con$padj),]                        # Removing genes with NA adj p-values                  

Visualizing results